A striped supersolid phase and the search for deconfmed quantum criticality in 

hard-core bosons on the triangular lattice 
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Using large-scale quantum Monte Carlo simulations we study bosons hopping on a triangular 
lattice with nearest (V) and next-nearest (V) neighbor repulsive interactions. In the limit where 
V — but V' is large, we find an example of an unusual period-three striped supersolid state that is 
stable at 1/2-filling. We discuss the relationship of this state to others on the rich ground-state phase 
diagram, which include a previously-discovered supersolid, a uniform superfiuid, as well as several 
Mott insulating phases. We study several superfiuid- and supersolid-to-Mott phase transitions, 
including one proposed by a recent phenomenological dual vortex field theory as a candidate for an 
exotic deconfmed quantum critical point. We find no examples of unconventional quantum criticality 
among any of the interesting phase transitions in the model. 
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I. INTRODUCTION 

The question of the effects of frustrated repulsive inter- 
actions on the Bosc Hubbard Hamiltonian is an intriguing 
one. Over the last several years, studies of the simplest 
hard-core boson Hamiltonian on two-dimensional (2D) 
frustrated lattices have revealed fascinating non-trivial 
quantum phases X The two most exciting recent examples 
are a stable supersolid phase in triangular lattice bosons 
at 1/2-filling )2iiLi and a valence-bond solid (VBS) phase 
at 2/3-filling on the kagome lattice^&i The VBS order 
corresponds to a short-range resonance of bosons around 
a lattice bond or plaquette, the supersolid phase to a full 
derealization of bosons across the lattice, resulting in 
off-diagonal long-range order at zero-temperature. 

The desire to extend these studies to other models in 
this general class is immediate. There is particular in- 
terest in the community to find other examples of non- 
trivial bulk quantum phases, and possibly to relate exist- 
ing ones to more realistic models of cold atoms trapped in 
optical lattices. The inclusion of further-neighbor inter- 
actions in the Bose Hubbard Hamiltonian is extremely 
important for efforts^ that attempt to use it to model 
a dilute gas of atoms with dipolar interactions (such as 
52 Cr )2 in a periodic potential. In addition, perturbations 
of the simplest hard-core boson Hamiltonian offer the 
opportunity to make connection with continuum quan- 
tum field theories studying the phenomenology of exotic 
quantum phase transitions. In this work, we enlist the 
results of a recent dual vortex theory (DVT)i£ to help 
perform a "guided search" for examples of exotic (non- 
Landau-Ginzburg- Wilson (LGW)) deconfmed quantum 
critical phenomenal 

We consider a model of bosons on the triangular lattice 
with nearest-neighbor (nn) hopping and both nearest and 
next-nearest neighbor (nnn) repulsion, 
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where b J (b { ) creates (destroys) a particle on site i, 
Tii = b\b i ,t> and V > denote the nn hopping and re- 
pulsion, V' > is a nnn repulsion, and \i is the chemical 



potential. In the usual way (Sf 



■1/2, S+ = b] etc.) 



the model can be mapped to the spin- 1/2 XXZ model on 
the triangular lattice, 
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where J± is the in-plane exchange, J z and J' z > denote 
the diagonal nn and nnn exchange, respectively, and h 
is an external magnetic field. The model parameters are 
related by Jj_ = 2t, J z = V, J' z = V, and h = n - W - 

w. 

This paper is organized as follows: In Section |TT] we 
review what is known about this class of Hamiltonian, 
including past work and general theoretical expectations. 
In Section IIIII we discuss recent progress and implemen- 
tations of quantum Monte Carlo (QMC) algorithms for 
the t-V-V'-fi model (Eq. {TJ)) taking special note of ef- 
forts in alleviating some common ergodicity problems. 
QMC results are presented in two sections: first, results 
at 1/2-filling (h — 0) for the t-V-V model, with a focus 
on the full phase diagram and in particular supersolidity 
in the nnn case (V = 0) (Section IIV[) . In this section 
we also study several superfiuid- and supersolid-to-Mott 
phase transitions with the aim of searching for examples 
of interesting quantum phase transitions. Second, results 
are presented in Section [V] for the V = 0, t-V'-fi model 
away from 1/2-filling. There we are able to access an ex- 
ample of a superfluid-to-Mott insulator quantum phase 
transition, which was suggested by a recent dual vortex 
study^ to be a candidate for exotic non-LGW type be- 
havior. In our case we find that this phase transition is 
strongly first-order. 
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FIG. 1: (color online) A finite-size 6x6 (N = 36) triangu- 
lar lattice (with periodic boundaries) with one of three nnn- 
interacting sublattices illustrated. Our conventions for prim- 
itive vectors are labeled ai and a.2- 



II. PAST WORK AND GENERAL 
CONSIDERATIONS 

Various limiting cases of the t-V-V'-fx model are well 
understood. When t dominates, the presence of a 
uniform superfiuid phase has been tacitly established. 
The simple nn Hamiltonian (with V' = 0) has re- 
cently been studied by several groups using numerical 
techniques ^dull where the most interesting feature of 
the phase diagram is a phase transition out of the uni- 
form superfiuid into a supersolid phase (with coexisting 
superfiuid and diagonal long-range order) at V/t = 9. 
This supersolid phase can be viewed as arising from an 
order-by-disorder mechanism, where an extensive degen- 
eracy of the classical ground states of the interacting 
part of the Hamiltonian is lifted by quantum fluctua- 
tions (hopping) J£ The extensively-degenerate manifold 
consists of spin configurations that allows for each trian- 
gular lattice plaquette to be minimally frustrated (one 
spin up and two spins down, or vice versa). In the su- 
persolid phase, these configurations obtain diagonal long- 
range order when perturbed by t. Alternatively, the su- 
persolid state may be viewed as resulting from doping 
the 1/3 or 2/3 (boson density) filled Mott phases that 
occur at > with additional superfiuid bosons. 

The purely nnn model, V = 0, is expected to share 
some analogies with the above nn model. In the t = 
classical limit the three nnn interacting sublattices are 
entirely decoupled, leading each nnn sublattice to re- 
alize the frustrated classical groundstate of the equiv- 
alent nn triangular Ising model. An interesting ques- 
tion that immediately arises is whether the inclusion of 
a nn hopping t will promote a similar order-by-disorder 
effect in the quasi-long-range (1/^/r) diagonal density- 
density correlations. Such a question is non-trivial, since 
although present in the aforementioned nn triangular lat- 
tice model, a supersolid phase is absent in the analogous 
kagome lattice model where 1 / \Jr correlations do not ex- 
ist in the frustrated classical manifold. 6 ' 7 We address this 
question using QMC simulations in Section HVl of this pa- 
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FIG. 2: (color online) The three sets of ordering wavevectors 
for large repulsive interactions discussed throughout the text. 



per. 

To gain insight into the types of order which may de- 
velop when large repulsive interactions (V and V) are 
present, we appeal to a simple calculation of the order- 
ing wavevectors for the classical nn and nnn Ising model. 
The "star" of k-points that minimize the Fourier trans- 
form of the classical interaction matrix are displayed in 
Fig. [2] and have the following values: for V' /V = 0, 

k 01 = (27r/3,-27r/V^) 
k a2 = (4tt/3,0), 

(and negatives, where the third can be constructed by 
adding a reciprocal lattice vector) which correspond to 
the corners of the first Brillouin zone; for V/V = 

k bl = (2tt/3, -2tt/3V3) 
k b2 = (0,4tt/3\/3) 
k b3 = (2^/3,2tt/3V3); 

and for V = V: 

k Cl = (tt, -7r/V3) 
k C2 - (0,2tt/V3) 

k C3 = (7T,7r/\/3). 

The k c are bisectors of the boundary of the first Bril- 
louin zone, and k b = 3/2k c . In the QMC, these ordering 
wavevectors are employed to construct diagonal order pa- 
rameters using the density structure factor, as discussed 
in the next section. 

In the presence of a moderate particle-hole symmetry- 
breaking chemical potential (magnetic field), the classi- 
cal nn model on the triangular lattice realizes a \/3 x \/3 
(boson density 2/3 or 1/3) pattern with true long-range 
ordered correlations at wavevectors k a J£ The same can 
be expected for the three nnn sublattices if not coupled 
by a hopping (i = 0). If one combines these three non- 
interacting nnn sublattices to form the original lattice 
illustrated in Fig. [fl two possible global patterns result, 
with ordering wavevectors k b . These are the "striped" 



3 



and "bubble" phases (illustrated in Fig.HO]). studied pre- 
viously using numerical techniques by Metcalf^ Met- 
calf demonstrated that both phases are stable, energet- 
ically equivalent ground-state configurations of the clas- 
sical nnn model in the presence of a magnetic field. We 
therefore can ask whether either of these ordered (or in- 
sulating) patterns survive in the system upon inclusion of 
the nn hopping t. If so, then the possibility of accessing 
the associated superfluid-insulating quantum phase tran- 
sition (say upon tuning t) is of exceptional interest. This 
phase transition was studied by Burkov and Balcnts 10 
using a phenomenological field theory of vortices inter- 
acting on the dual (honeycomb) lattice, with the goal 
of understanding universal phenomena in the vicinity of 
various superfluid-Mott quantum critical points in mi- 
croscopic models with the same spatial symmetries. The 
simplest DVT analysis for the case of 1/3 (or 2/3) bo- 
son filling contains three unique Mott insulating states. 
Two of these Mott phases are the striped and bubble 
states, and are conjectured to be connected by a decon- 
fined quantum critical point (DQCP) to the neighbor- 
ing superfluid state. A DQCP is the recently proposed 
bon ton of unconventional quantum criticality, not based 
on the standard LGW concept of the order parameter fii 
Rather, such DQCPs are described in terms of emergent 
degrees of freedom not related to the microscopic sym- 
metries of the adjoining phases, and emergent topological 
order. However, since no clear examples of a DQCP have 
been identified in microscopic models or real materials, 
the discovery of the relevant superfluid-insulating quan- 
tum phase transition in our model would offer the unique 
opportunity to search for this exotic phenomenon guided 
by predictions of the same phenomenological field theory 
from which it was originally proposed, — In Section fVl we 
report QMC evidence for the existence of the striped and 
bubble phases in the t — V — fi model and study in detail 
the relevant superfluid-insulating phase transition. 



III. NOTES ON QMC TECHNIQUES 

In an effort to elucidate the phase diagram of the 
model and study the physical questions mentioned above, 
we employ extensive numerical simulations in this work, 
using the Stochastic Series Expansion (SSE) quantum 
Monte Carlo method pioneered by Sandvik i 17 ' 18 The tri- 
angular lattice with next-nearest neighbor interactions is 
illustrated in Fig. [TJ Here, the nnn diagonal interactions 
form three decoupled triangular sublattices, which are 
coupled by the nn hopping term. On frustrated lattices 
such as this, the QMC is free from the sign problem only 
for positive (ferromagnetic) t. The sign of V or V' is not 
relevant in this regard, however as observed in studies on 
other frustrated lattices, simulation efficiency is severely 
affected by large repulsive interactions. Such problems 
arise due to the difficulty of tunneling between the classi- 
cally degenerate manifold of minimally-frustrated states; 
in classical simulations of such systems, non-local clus- 



ter or loop moves are often employed in a Monte Carlo. 
However for this type of finite-temperature QMC, the 
analogous cluster moves are highly non-trivial to imple- 
ment. 

One significant algorithmic advance, which combats 
the tendency of the SSE directed loop algorithm to freeze 
(for large V or V') is to employ six-leg vertices rep- 
resenting triangular lattice plaquettes instead of 4-leg 
lattice bonds in the formulation of the directed loop 
algorithm,!^! 2 ! In the traditional SSE formulation,! 7 .!! 8 , 
the Hamiltonian is decomposed into a list of two-site 
bond operators, which are sampled using diagonal and 
off-diagonal updates in a Monte Carlo procedure. The 
off-diagonal update typically employed is an operator 
loop algorithm called the "directed loop" update: closed 
loops are constructed in the d + 1 simulation cell ac- 
cording to transition weights determined by the partition 
function. After being built, bosons on these closed loops 
are added or removed with probability 1. 

A difficulty in this traditional bond algorithm occurs 
for certain parameters, particularly in the large V or 

V regime of the phase diagram (the frustrated near- 
degenerate manifold with small hopping). As described 
in Ref. 1 2 11 when the loop under construction encounters 
a frustrated bond with some diagonal operator, the non- 
zero transition weight associated with adding or remov- 
ing a single boson (resulting in the sampling of another 
diagonal operator) is suppressed due to the repulsive in- 
teractions. 

This difficulty can be overcome by a removal of the 
"local" energy barrier directly associated with this sup- 
pression of transition weight. To do so, one may consider 
a decomposition of the Hamiltonian into a list of three- 
site triangular plaquettes. The plaquette algorithm can 
be expected to increase the efficiency of the directed- 
loop updates; movement between different minimally- 
frustrated plaquettes are no longer suppressed because 
both the initial and final configurations of the plaque- 
tte have the same (or similar) transfer weights. The 
result is a removal of an energy barrier directly asso- 
ciated with traversing local configurations in the near- 
degenerate frustrated manifold for states located in the 
large V or V' regime of the phase diagram, and a signif- 
icant improvement in efficiency for the algorithm. 21 

We employ the triangular plaquette version of the SSE 
algorithm for both nn and nnn Hamiltonian terms. As 
discussed below, we nonetheless observe "freezing" of the 
QMC dynamics for V'/t > 5, due to the persistence 
of "global" energy barriers associated with the near- 
degeneracy of the entire minimally-frustrated manifold 
and not alleviated by the above algorithmic advances. 
This freezing problem could be further reduced for larger 

V or V' by employing more advanced simulations tech- 
niques, such as parallel tempering^ however, we find 
that the essential physics is captured for the parameter 
range accessible in this simulation scheme, so we leave 
such improvements for future studies. 

The standard estimator for the superfluid density of 
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the model is the stiffness^ defined via the free energy 
response to a twist <f> in the boundary conditions and 
measured with winding number fluctuations^ in each of 
the lattice translation vector directions, 

a-Sf. (3) 

where (3 — t/T is the inverse temperature. Typically 
the stiffness is averaged over both ai and &2 directions 
(see Fig. [1]) so that we define p s = (p* 1 + p* 2 )/2, unless 
measured in a phase which breaks rotational symmetry 
(we discuss this more in the next section). 

The presence of a non-zero spin stiffness in the ther- 
modynamic limit indicates off-diagonal long-range or- 
der, manifest as either a uniform superffuid or supersolid 
phase. The standard indicator for the presence of diago- 
nal long-range order is the density structure factor, 

= ^E eJk ' (rs ~ rj) ((^> - (4) 

easily measured in the rii basis, which is diagonal in our 
SSE implementation. The order parameter for diagonal 
correlations is constructed from the structure factor by 
noting the non-zero ordering wavevectors k Qi through the 
definition 

«> = % 1 - (5) 

where N is the finite number of lattice points. Any stable 
phase having non-zero values for both p s and (m^.) in 
the thermodynamic limit can be classified as a supersolid, 
while we refer to the phase with only (m^ . ) ^ as a Mott 
insulator. 



IV. HALF-FILLING 

A. Next-nearest neighbor striped supersolid 

We begin by studying the V — Hamiltonian at 1/2- 
filling, or /i = 3V' (corresponding to zero magnetic field 
in the XXZ model). As described in Section [Til m the 
classical limit (t — 0) the existence of the nnn V' inter- 
action alone produces three decoupled triangular lattices, 
each of which realizes the usual ground state described by 
the set of wavevectors . It is most easily described as 
a degenerate manifold of minimally-frustrated triangular 
plaquettes, each with only one or two bosons. It is ex- 
pected that the addition of a nn hopping (t) will break the 
classical degeneracy, although whether this perturbation 
produces an ordering effect on the diagonal components 
of the bosons is not obvious. 

We can address this physical question with QMC sim- 
ulations, however, pragmatically, the approach of per- 
turbing the classical V manifold with t is difficult due to 
the aforementioned severe algorithmic freezing effects in 




FIG. 3: (color online) Superfluid-supersolid phase transi- 
tion. The top panel illustrates the superfluid density (p s — 
(ps 1 + Ps 2 )/2), and the bottom the solid order parameter, as 
a function of the nnn repulsion. Data for both main plots 
has been extrapolated to the thermodynamic limit. The in- 
set shows a detail of the superfluid density in the supersolid 
phase. The values pi and pjr are the superfluid density paral- 
lel and perpendicular to the stripe direction, discussed in the 
text. 



this regime. We thus employ the usual complementary 
approach, by applying a strengthening diagonal V inter- 
action to the dominant hopping regime (we use a fixed 
t = 0.5 corresponding to J± = 1 throughout this paper). 
This allows us to study the effects of increasing frustra- 
tion while simultaneously evaluating the simulation per- 
formance in a systematic way. Similarly, ground-state 
estimators are obtained by the usual method of cooling 
the simulation cell (by lowering T systematically) until 
convergence is reached to within error bars. Since al- 
gorithmic freezing strongly affects error bars, a delicate 
balance must be struck in order to find a temperature 
that faithfully represents the ground-state physics, while 
at the same time retaining ergodicity on reasonable sim- 
ulation time-scales. 

Results of the simulations are illustrated in Fig. [3] 
As one increases the frustrating V', we observe a phase 
transition out of the usual "uniform" superfluid phase at 
dominant t. This phase transition is indicated by a drop 
in the p s curve in Fig. [3] at V' ~ 3, while p s remains finite 
(in the thermodynamic limit) on both sides of the tran- 
sition. Note that p s is measured in units of t, which has 
been set to 0.5 throughout this work. For V > 3, long- 
range order develops in the diagonal density correlations, 
as indicated by very sharp Bragg peaks in the structure 
factor S'(k) at ordering wavevectors which correspond 
to the expectations from arguments involving the sym- 
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FIG. 4: (color online) Examples of finite-size scaling for quan- 
tities that appear in the main plot of Fig. [3] for two values of 
V. 



metry of the interactions in the nnn Ising model (Sec- 
tion Hi). The diagonal long range order is classified using 
the set of order parameters (m^.) (Eq. (|5)). The large V 
phase exhibits concurrent diagonal long range order (as 
indicated by the survival of the order parameter (m%} in 
the thermodynamic limit) and off-diagonal long range or- 
der (indicated by a non-zero superfiuid density) signaling 
that the entered phase is a thermodynamically stable and 
homogeneous supersolid at 1/2-filling. The survival of p s 
and (ml) are rigorously confirmed using finite size scal- 
ing, illustrated in in Fig. [U We observe that p s quickly 
converges to its thermodynamic limit value, while (m%) 
requires the simulation of relatively large lattices in the 
supersolid phase - this apparent non-monotonicity in the 
scaling is discussed further below. In the uniform super- 
fluid phase, (ml) extrapolates to zero as expected for a 
phase without diagonal long range order. Deeper in the 
supersolid phase, for V' > 5, we observe that simulation 
efficiency is quickly destroyed, manifest as observations 
of non-equilibrium freezing in the measured observables. 

The ordering wavevectors at are consistent with 
a striped or bubble phase as observed in the classical 
nnn Ising model (see Section |TT), and could be expected 
in a canonical simulation fixed at 1/3 or 2/3- filling. 16 
Our simulations are however grand-canonical and the ob- 
served supersolid is firmly at 1/2-filling, as indicated by 
studying histograms of the particle density. The super- 
solid is homogeneous and shows no obvious evidence of 
phase separation. 

In a phase such as a supersolid, with hopping bosons 
delocalized throughout the lattice, one cannot distin- 
guish the particular diagonal sublattice ordering in QMC 
simulations by looking directly at density "snapshots" . 
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S(k b ,)S(k bj ) 

FIG. 5: (color online) The three nnn striped diagonal order 
parameters (5(kb ; )) plotted as a function of Monte Carlo time 
(top) for a 30 x 30 lattice. Each data point is a bin average 
of 10 4 QMC steps. A histo gram of the three possible product 
combinations (bottom) for single bin of 10 4 QMC steps indi- 
cates the presence of striped correlations, as discussed in the 
text. The stiffness values calculated for this same simulation 
were p? 1 = 0.398(8) and pf 2 = 0.384(8). 



Any obvious hallmarks of ordering patterns tend to be 
washed out by the hopping. We can however learn more 
about the sublattice ordering by studying the symmetry 
of the ordering wavevectors. Specifically, in the thermo- 
dynamic limit, the striped phase (see also Fig. flO) is char- 
acterized by the existence of only one nonzero value for 
(m?), whereas the bubble phase has all three components 
nonzero and equal in magnitude. Monte Carlo data for 
the three independent 5(kh) is illustrated in Fig. [5] for 
a large system size in the supersolid phase. It is clear 
that the three order parameters "trade off" between a 
large and small value. If one remembers that S(k&(£)) 
retains some non-zero magnitude even in a disordered 
phase due to finite size effects, it is apparent that the 
behavior in Fig. [5] suggests a stripe-like symmetry. Fur- 
ther evidence can be obtained by looking at histograms 
of the products of the order parameters: S'(k(, 1 )S , (kh 2 ), 
S , (kb 2 )5(kb 3 ), and S l (k(, 1 )S , (k(, 3 ). In a finite-size remnant 
of the striped phase, one expects one of these products 
to be small and the other two large. The bottom plot in 
Fig. [5] is clearly consistent with this expectation. With 
this in mind, we note that one must be careful to include 
only the largest 5(k&) estimator in the construction of the 
striped order parameter (for example in Fig. [3) . We ob- 
serve that below a system size of 12 x 12, no single element 
of S^k;,) seems to be preferred. The non-monotonicity in 
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the finite-size scaling of (mT) (for V' — 3.6) in Fig. [4]ap- 
parently occurs at the point where the finite simulation 
cell becomes large enough to see the length scale asso- 
ciated with stripe development. This interesting feature 
should be investigated further. 

The above evidence points to the naive picture of the 
supersolid phase as a doped version of the 1/3 or 2/3 
period-three striped Mott phase. However some of the 
behavior of the supersolid is unusual in this regard. In 
particular, as alluded to above (see Fig. [5]), p s does not 
show strong anisotropy or correlation with stripe "direc- 
tion" (i.e. the largest value of S(k bl ), S(k b2 ) or S(k b3 )) 
such as one might expect for a system with broken ro- 
tational symmetry. To study this further, we observe p s 
along both primitive lattice directions during the QMC 
run. In a strongly developed striped phase with large 
S(k bl ), one might expect the parallel-running superfluid 
density p^ 2 to be largest, while in a S(k b2 ) striped phase, 
Ps 1 might dominate. We ran simulations that, for each 
QMC configuration sampled, binned into a new esti- 
mator p" s or pjr. If stripes running parallel to ai were 
observed (i.e. Sik b2 ) is largest), p" s = p* 1 and pjr = pf 2 ; 
vice versa for stripes running parallel to a.2- If stripes 
were observed to run parallel to the third symmetry di- 
rection, i.e. S^kb,) is largest, no data for p" s or pjr was 
binned, since one might expect the superfluid density to 
be isotropic in this case. 

We explore the possibility for p s anisotropy in the 
striped supersolid phase. Data for a region of the su- 
persolid phase adjacent to the transition is illustrated in 
the inset to Fig. [31 taken for a 30 x 30 system at (3 = 12.5 
This measurement indicates that statistically significant 
anisotropies are present, with deviations of a few per- 
cent, for example at V' = 3.4, where p" s — 0.484(2) 
and pjr = 0.467(3). However, it is clear that large 
anisotropies like those reported in different striped su- 
persolid states in Bose Hubbard models on the square 
lattice^ are not evident in our nnn triangular striped 
supersolid. This suggests that the usual interpretation 
of the supersolid as a doped Mott insulator may not be 
applicable here. We revisit this discussion after looking 
at the period-three striped Mott phase in Section IVl 



B. V-V' phase diagram 

Now that we have firmly established the presence of 
a nnn striped supersolid phase at large V that exists 
complementary to the previously discovered nn super- 
solid, it is natural to extend our exploration of the phase 
diagram by perturbing into the V, V > region. An 
interesting Mott state is uncovered, which exists at 1/2- 
filling for moderate to large V, V > 0. The phase is 
a period- two striped state (see Fig. [5]) , with boson den- 
sity modulated at one of the ordering wavevectors k c = 
(tt, -71-/73), (0,271-/^(3) or (7r,7r/73). These wavevec- 
tors are the same as those expected from the minimiza- 
tion of the ground state energy for the classical (t = 0) 




FIG. 6: (color online) Schematic ground-state phase diagram 
for the zero-field t-V-V' model. Blue (dashed) lines represent 
transitions between the superfluid and supersolid states. The 
red (thick) line is the phase boundary between these phases 
and the period-two Mott insulator, which is strongly first or- 
der. Dotted lines labeled I, II and III illustrate the position 
of data sweeps in the next three figures. The inset shows 
one of the possible period-two striped Mott phases, found by 
plotting the boson density for an instantaneous QMC config- 
uration with V = 5.0 and V = 0.8. 



limit of the model in this parameter regime, as discussed 
m Section HH 

The phase boundaries between this Mott state and the 
superfluid and supersolid phases in Fig. [6] are interest- 
ing examples of quantum phase transitions that we can 
characterize using the QMC simulation. Data along three 
cuts in the phase diagram in Fig. [S] are discussed in detail 
below. 

First, we examine cut I by starting the simulation in 
the nn supersolid phase (at a value of V = 5) and "per- 
turbing" away from the x-axis by increasing V' (Fig. [7]) . 
It is evident upon increasing V' that the magnitude of 
the nn diagonal order parameter ((fn^)) decreases until 
terminating in a phase transition, upon which p s falls 
abruptly to zero and the period-two Mott order param- 
eter ((m^.)) jumps almost to its maximum value. The 
sharp discontinuity at the phase boundary is fairly typi- 
cal of a strong first-order transition. 

A qualitatively similar behavior occurs along cut II as 
we perturb the nnn striped supersolid at V = 3.4 by the 
interaction V. In this case, a much smaller lattice was 
employed due to the observation of severe algorithmic 
freezing in the simulation for moderate to large values of 
V and V. Nonetheless, the abrupt jumps in the observ- 
ables in Fig. [5] are suggestive of a finite-size remnant of a 
strong first-order phase transition, even for a system size 
as small as 6 x 6. 

Finally, we examine the remaining region of the phase 
boundary by beginning the simulation in the uniform su- 
perfluid phase (V = V = 0) and perturbing by increas- 
ing V = V. We again encounter an apparent direct first 
order phase transition between the superfluid and period- 
two Mott state, with no other obvious intervening phases 



FIG. 7: (color online) Superfluid density, nn supersolid order 
parameter, and period-two Mott order parameter along cut I 
in the phase diagram. Data was taken at V = 5 for a L = 12 
lattice at /3 — 12.5. 



FIG. 9: (color online) Superfluid density and period-two Mott 
order parameter along cut III in the phase diagram. Data was 
taken for a L — 18 lattice at /3 = 12.5. 




FIG. 8: (color online) Superfluid density and period-two Mott 
order parameter along cut II in the phase diagram. Data was 
taken for a L = 6 lattice at V' — 3.4 and = 10. 



occurring (Fig. [5]). 

The superfluid- and supersolid-to-Mott phase bound- 
aries in Fig. [5] appear to be simple, strongly first-order 
phase transitions. This result is perhaps expected in a 
typical LGW paradigm since very different symmetries 
are broken on each side of the phase boundary. In con- 
trast, we turn now to an example of a superfluid-to-Mott 
quantum phase transition where the theoretical possibil- 
ity exists for exotic non-LGW behavior. 



V. AWAY FROM 1/2-FILLING, fx > 3 AND V = 




•••••••• 




FIG. 10: (color online) QMC "snapshots" of instantaneous 
boson density (spin) arrangements stabilized on a 18 x 18 
simulation cell for (3 = 12.5, V = 2.75 and y, = -4.5. The 
defects seen in these pictures are not static. Configurations 
shown have density 2/3, and are referred to in the text as the 
"bubble" (top) phase and "striped" (bottom) phase. 



Consider the model Hamiltonian Eq. ([T]) with the nn 
repulsion set to zero (V — 0), and V large enough to 
be inside the nnn supersolid phase. The application of 
a chemical potential to tune the supersolid phase away 
from 1/2- filling reveals interesting features. The phase 
diagram is observed to have the usual "lobe" structure, 
characteristic of models with frustrated boson repulsion. 
Explicit illustrations of this phase diagram can be found 
in Refs. 0, and [l5|. Of particular interest are the two 



lobes of Mott insulator at fillings 1/3 and 2/3, which in 
the present case are separated by the striped supersolid 
state. These Mott insulating regions are observed to have 
the same ordering wavevectors as the supersolid (kj). 
Further, the bosons in these Mott phases are strongly 
localized, allowing one to observe the sublattice density 
directly from "snapshots" of the diagonal QMC simula- 
tion basis, something that is not possible in the supersolid 
phase. 
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FIG. 11: (color online) The superfluid-Mott transition at 2/3 
filling, for a L — 18 system. The simulation parameters were 
13 = 12.5, and h = 3.750 (ji = h - W). This value of h was 
found to very closely match the filling of the superfluid right 
at the transition to 2/3. The main plot illustrates the abrupt 
discontinuity in p s characteristic of a first-order transition. 
The inset shows hysteresis typical of a first-order energy level 
crossing. 



Interestingly, in the simulations the Mott insulating 
phases appear able to manifest both expected (period- 
three striped or bubble) sublattice ordering patterns, de- 
pending on equilibration details, and what starting con- 
figurations were used in the QMC. Such behavior oc- 
curs due to an algorithmic non-ergodicity that may cause 
"sticking" , i.e. a suppression of free sampling in these two 
low-lying states. Examples of the states are illustrated 
in Fig. 1 1 01 where two independent simulations were run 
for parameters inside of the 2/3 Mott lobe starting from 
different random basis configurations. One simulations 
settles firmly in the striped state, the other in the bub- 
ble state. After fairly long run-times (10 8 QMC steps), 
we do not observe any "mixing" between states, i.e. one 
simulation remains in the striped state and the other re- 
mains in the bubble state. Although boson hopping pro- 
cesses are not entirely suppressed (indeed the "defects" 
illustrated in Fig. [TU] are not fixed), large energy barriers 
exist that prevent these two examples of sublattice order- 
ing from mixing. The preference of the model towards 
either state may be delicate, or vary with position on the 
phase diagram. Nonetheless we observe that the average 
energy of the striped state E — —2.72306(4) (measured 
in units of t = 0.5) is slightly lower than the energy of 
the bubble state E — —2.72280(4) for the specific param- 
eters discussed in Fig. [TU] It is interesting to note that 
a preference of the system towards the 2/3 striped state 
is consistent with the picture of this phase being closely 
connected to the adjacent nnn supersolid phase, i.e. upon 
crossing the Mott to supersolid phase boundary, the bro- 
ken symmetry of the diagonal order parameter stays the 
same. 

As discussed in Section QTJ a recent DVT— predicts 
that two specific 2/3 Mott states are connected to a 



neighboring superfluid (at the same filling) by an exotic 
DQCP. These two Mott states are the striped and bubble 
states, both of which are found in our QMC simulations 
as described above. We can easily access the superfluid- 
Mott phase boundary by varying model parameters (h or 
V') in the QMC. Thus, with our simulation we are in the 
unique position of possibly making a direct connection 
between the behavior of a simple microscopic model and 
an exotic phenomenon contained in a continuum DVT. 

The first step in the simulations is to study the order 
of the phase transition. The "conventional" situation in 
a case like this where the phases on either side of the 
phase transition have unrelated broken symmetries is a 
first-order transition; a continuous quantum critical point 
would suggest a breakdown of the conventional Landau 
picture and possible relevance to a more exotic theory 
like deconfined quantum criticalityjii 

Using the QMC simulation, we mapped out the phase 
diagram roughly near the "tip" of the Mott lobe, in at- 
tempt to determine a set of parameters which closely 
matches the filling of the neighboring superfluid phase to 
2/3 (which would be required to make a comparison to 
the predictions of the fixed-filling DVT). We then per- 
formed more detailed simulations, keeping the chemical 
potential fixed and sweeping the parameter V'. The re- 
sults are illustrated in Fig. [TTJ These simulations, even 
though carried out on a relative small lattice (18 x 18), 
display a discontinuity in the value of p s and hysteresis 
effects in other observables - behavior that is prototypical 
of a first-order phase transition. In fact, the superfluid- 
Mott phase transition appears "strongly" first-order, e.g. 
the discontinuity in p s is clear and obvious. This result 
precludes an interpretation of the transition as a quan- 
tum critical point that should share universal properties 
with a DQCP, and instead suggests the much more con- 
ventional LGW paradigm as the correct descriptor of the 
transition. 



VI. DISCUSSION 

In this paper we have presented QMC simulation re- 
sults for the ground state phase diagram of the Bose 
Hubbard model with nn hopping and nnn repulsive in- 
teractions. At 1/2-filling, the phase diagram contains a 
superfluid phase (for dominant hopping (t)), two super- 
solid phases (for dominant V and V respectively), and 
a period-two striped Mott insulating state for large V 
and V. The superfluid- and supersolid-to-Mott phase 
boundaries are observed to be strongly first-order. 

The large- V supersolid phase is characterized by the 
homogeneous coexistence of diagonal (solid) and off- 
diagonal (superfluid) long range order in the system. 
Density correlations in this phase are observed to have 
the symmetry of a period-three striped state, similar to 
that observed in an adjacent 1/3-filled striped Mott in- 
sulator, although the density of the supersolid is mea- 
sured to be firmly at 1/2-filling (for h = 0). This fact 
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leads to the immediate interpretation of the supersolid 
state as a 1/3 Mott insulator which is doped with in- 
terstitial bosons that become delocalized, forming the 
"superfiuid" component. This picture is clearly an over- 
simplification, however, since as a result one would ex- 
pect the boson superfiow to occur in one-dimensional 
channels between the strongly-localized density stripes. 
Indeed the resulting anisotropy in the superfiuid density 
is measured, although it is extremely small (w 5%), espe- 
cially when compared to other striped supersolid states 
observed on square-lattice Bose Hubbard models doped 
away from 1/2-filling^ This issue might be answered by 
looking at correlations further away from the superfluid- 
supersolid quantum critical point, or examining other es- 
timators that would be expected to harbor a comple- 
mentary anisotropy in such a state. Since it is also possi- 
ble that the doped-Mott insulator interpretation simply 
ceases to be valid at densities this far from 1/3 or 2/3, 
additional measurements of the p s anisotropy away from 
1/2 filling would help shed light on the issue. The reso- 
lution of such paradoxes in a simple model like this could 
have broad implications for studies of nanoscale density 
inhomogeneities in other bulk quantum phases, for exam- 
ple in "stripe" mechanisms of high-T c superconductivity 
in the cuprates. 

The 1 /3-filled striped Mott insulating state exists on 
the V — phase diagram for large V, in the presence of 
a chemical potential that breaks particle-hole symmetry. 
In addition to being adjacent to the striped supersolid 
phase, this phase borders on a uniform superfiuid phase 
that dominates at larger hopping t. This phase boundary 
offers a unique opportunity to search for an example of 
exotic non-LGW quantum criticality suggested by a re- 
cent phenomenological DVT. 10 However, in this model, 
the relevant superfluid-Mott phase transition is observed 
to be strongly first-order. Thus an example of a highly 
sought deconfined quantum critical points is not present. 
Although it is true that the DVT is a phenomenological 
theory and does not predict the order of a phase tran- 
sition for any specific microscopic model, it is perhaps 



unfortunate that the expected universal properties of the 
vortex theory are not realized in a simple model such as 
this. 

On the other hand, the result that the relevant 
superfluid-Mott phase transition is strongly first order is 
an important contrast to other phase transitions previ- 
ously studied as candidates for DQCPs. In particular, 
simulations of superfiuid- VBS phase transitions 7,25 which 
at first appear to be continuous (for example because of 
the lack of obvious discontinuities in the data) may in- 
deed turn out to be very "weakly" first-order upon closer 
scrutiny of more subtle indicators, such as values of crit- 
ical exponents or interpretations of finite-size behavior 
on extremely large lattices. In contrast to these previ- 
ous simulations that involved transitions into the VBS 
phase, our present candidate transition only requires the 
presence of a simple density- modulated Mott insulator. 
Comparing the numerical signatures between these cases 
may help develop a more systematic method for iden- 
tifying weakly-first order phase transitions in quantum 
simulations, which could contribute to the resolution of 
fundamental controversies 26 surrounding the very exis- 
tence of DQCPs in these and related systems. 
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